{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inverse Transform Sampling\n",
    "\n",
    "This notebook will conceptualize how inverse transform sampling works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy.random as ra\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a spectrum which follows an `almost` bell-curve type distribution (anyway, the specific type of distribution is not important here). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 2, 3, 4, 5, 6], [2000, 4040, 6500, 6000, 4020, 2070]]"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spectrum = [[1, 2, 3, 4, 5, 6],[2000, 4040, 6500, 6000, 4020, 2070]]\n",
    "energies = np.array(spectrum[0])\n",
    "fluxes = np.array(spectrum[1])\n",
    "spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, first we compute probabilities of flux. Afterwards, we compute the cumulative probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.08120179,  0.2452294 ,  0.5091352 ,  0.75274056,  0.91595615,  1.        ])"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob = fluxes/float(sum(fluxes))\n",
    "cum_prob = np.cumsum(prob)\n",
    "cum_prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We draw ten thousand numbers from uniform random distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.49834338,  0.31993222,  0.35882619,  0.15837646,  0.22595417,\n",
       "        0.85575223,  0.85203039,  0.78380252,  0.04170078])"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10000\n",
    "R = ra.uniform(0, 1, N)\n",
    "R[1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assign energies to events corresponding to the random number drawn.\n",
    "\n",
    "_Note: The command below finds bin interval using a single command. I am not sure though that it's very readble. Would\n",
    "we want to split that in multiple lines and maybe use explicit loops to make it more readable? Or is it fine as it is?\n",
    "Comments?_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 3, 3, 2, 2, 5, 5, 5, 1]"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gen_energies = [int(energies[np.argwhere(cum_prob == min(cum_prob[(cum_prob - r) > 0]))]) for r in R]\n",
    "gen_energies[1:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram energies to get shape approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 825, 1652, 2626, 2466, 1589,  842], dtype=int64)"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gen_energies = ((np.array(gen_energies) - 1) / 1).astype(int)\n",
    "times = np.arange(1, 6, 1)\n",
    "lc = np.bincount(gen_energies, minlength=len(times))\n",
    "lc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEPCAYAAAC3NDh4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0VNXXxvHvDoROQgnSiRQpiihgAREMJaEX6UGqoqiI\nYgcLYEHUH/gqRSnSBWkiIDWAhN5BRKmCgDSlhg4h2e8fGWLADEwgk5uyP2vNyszcMs8MYXbOPfee\nI6qKMcYYEx8fpwMYY4xJvqxIGGOMccuKhDHGGLesSBhjjHHLioQxxhi3rEgYY4xxy+tFQkTqiMgO\nEdklIm/Hs7yRiGwRkc0isk5Eqni6rTHGGO8Sb14nISI+wC6gJnAYWA+0VtUdcdbJoqoXXPfvB6ao\nahlPtjXGGONd3m5JPALsVtX9qhoJTAIax13hWoFwyQZEe7qtMcYY7/J2kSgI/BXn8UHXc9cRkSYi\nsh34CXg6IdsaY4zxnmTRca2qM1S1DNAE+NjpPMYYY2Kk9/L+DwFF4jwu5HouXqq6QkSKiUiuhGwr\nIjYAlTHGJJCqyq3W8XZLYj1QQkQCRSQD0BqYFXcFESke534FIIOqnvRk27hU1W6q9O7d2/EMyeFm\nn4N9FvZZ3PzmKa+2JFQ1SkReAsKIKUgjVXW7iHSJWazDgWYi0h64AlwEWt5sW2/mNcYYcz1vH25C\nVecDpW54blic+58Dn3u6rTHGmKSTLDquTeIJCgpyOkKyYJ/Dv+yz+Jd9Fgnn1YvpkoqIaGp4H8YY\nk1REBPWg49rrh5uMMe7dfffd7N+/3+kYJhULDAxk3759t729tSSMcZDrrzmnY5hUzN3vmKctCeuT\nMMYY45YVCWOMMW5ZkTDGGOOWFQljTJIpW7Ysy5YtS/LX7dSpE7169Ury102oiRMnUqdOHadjXMeK\nhDHmpoKCgsiVKxeRkZF3vK/ffvuNatWqJUKq1KlNmzbMnz/f6RjXsSJhjHFr//79rFixAh8fH2bN\ncjt0mvFAVFSU0xFuixUJY4xb48aNo3LlynTs2JExY8Zct2zu3Lncd999+Pn5UbhwYb744gsATpw4\nQcOGDcmZMye5c+fmiSeeiN2maNGi/PzzzwBcunSJDh06kCtXLu677z7+97//Ubhw4evWHTBgAA88\n8AA5c+YkNDSUK1euxC6fPXs25cuXJ2fOnDz++ONs3bo1dtnmzZupWLEi/v7+tG7dmkuXLt30fY4a\nNYp7772X3LlzU7duXQ4cOBC7zMfHh2HDhlGyZEly5crFSy+9lKBtv/76a0qWLEnJkiUBCAsLo3Tp\n0uTMmZOuXbsSFBTEqFGjABg7dixVq1aN3X7Hjh2EhISQO3duypQpw9SpU2/5+Sc6p0ciTKTRDNWY\nlCi5/+6WKFFChw4dqhs3blRfX1/9559/Ypflz59fV65cqaqqp0+f1s2bN6uqas+ePfWFF17QqKgo\nvXr1qq5YsSJ2m7vvvlsXL16sqqpvv/22BgUFaUREhB46dEjLlSunhQsXvm7dRx99VI8ePaqnTp3S\nMmXK6LBhw1RVddOmTXrXXXfp+vXrNTo6WseNG6d33323XrlyRa9cuaKBgYH61Vdf6dWrV3XatGnq\n6+ur77//frzvccaMGXrPPffozp07NSoqSvv27auPPfZY7HIR0YYNG+qZM2f0wIEDmidPHl2wYIHH\n24aEhOjp06f10qVLevz4cfXz89MZM2ZoVFSUfvXVV5ohQwYdOXKkqqqOGTNGq1atqqqq58+f18KF\nC+vYsWM1Ojpaf/nlFw0ICNDt27ff9PO/kbvfMdfzt/5+9WSl5H5L7v/RjHHnlr+7vXvH/De98da7\nt+fru1v3FpYvX64ZMmTQkydPqqpqmTJl9Msvv4xdHhgYqMOHD9czZ85ct12vXr20SZMm+scff/xn\nn3GLRLFixXThwoWxy7799tv/FImJEyfGPn7rrbf0hRdeUFXVF154QXv16nXdvkuVKqXLli3TZcuW\nacGCBa9b9thjj7ktEnXr1tVRo0bFPo6KitIsWbLogQMHVDXmi37VqlWxy1u2bKmfffaZx9uGh4fH\nLh83btx1RURVtXDhwvEWicmTJ2u1atWuW7dLly764Ycfqqr7z/9Gd1ok7HCTMclZnz7xlYiY5z1d\n3926tzBu3DhCQkLImTMnAKGhoYwdOzZ2+Q8//MCcOXMIDAykevXqrFmzBoC33nqL4sWLExISQokS\nJfjss8/i3f/hw4cpVKhQ7OO4h5quyZs3b+z9LFmycO7cOSCmr2TAgAHkypWLXLlykTNnTg4ePMjh\nw4c5fPgwBQteP9NxYGCg2/e5f/9+Xnnlldh95c6dGxHh0KF/5zi7WY5bbRv3PR4+fPg/7zPu8htz\nrVmz5rr3OHHiRP7++2/A/eef2GzsJmPMf1y6dIkpU6YQHR1N/vz5Abhy5QqnT59m69at3H///VSs\nWJEZM2YQFRXFoEGDaNmyJQcOHCBr1qz079+f/v37s23bNqpXr84jjzxC9erVr3uN/Pnzc/DgQUqX\nLg1w3bH8WylcuDDvvvsuPXv2/M+yZcuWXfclfW3fJUqUiHdfRYoU4b333iM0NNTj14+b41bbivw7\n8kX+/Pn/cwLAwYMH3e47KCiIBQsWxLvc3eef2KwlYYz5jx9//JH06dOzfft2tmzZwpYtW9i+fTtV\nq1Zl3LhxXL16lYkTJ3LmzBnSpUtH9uzZSZcuHQBz5sxhz549AGTPnp306dPHLourZcuW9OvXj9On\nT3Po0CGGDBnicb5nn32WoUOHsm7dOgDOnz/P3LlzOX/+PJUrVyZ9+vQMGjSIq1evMn369Nj14tOl\nSxc++eQTtm3bBkBERATTpk3zKMfzzz+foG3r16/Pb7/9xqxZs4iKimLw4MGxLYMbNWjQgF27dvHd\nd99x9epVIiMj2bBhAzt27CAyMtLt55/YrEgYY/5j3LhxPP300xQsWJC77ror9ta1a1cmTJgAwPjx\n4ylatCg5cuRg+PDhTJw4EYDdu3dTq1YtsmfPTpUqVejatWvstRFx/6ru1asXBQsWpGjRooSEhNCi\nRQsyZswYuzzuujeqWLEiI0aM4KWXXiJXrlyULFky9lCYr68v06dPZ/To0eTOnZupU6fSrFkzt/tq\n0qQJPXr0oHXr1uTIkYNy5cpdd63CjTniPk7ottfyvPnmmwQEBLBjxw4eeuih6973NdmyZSMsLIxJ\nkyZRoEABChQoQI8ePWLP8HL3+Sc2GwXWGAfZKLD/Gjp0KJMnT2bJkiVOR0kyqkqhQoWYOHHidacK\nJyYbBdYYkyIdPXqUVatWoars3LmTAQMG0LRpU6djeV1YWBgRERFcvnyZvn37AlCpUiWHU7lnHdfG\nGEdcuXKFLl26sG/fPnLkyEFoaCgvvPCC07G8bvXq1bRp04bIyEjuvfdeZs6cGe/hpuTCDjcZ4yA7\n3GS8zQ43GWOM8RorEsYYY9yyImGMMcYtKxLGGGPcsiJhjDHGLSsSxphE1a9fP5577rlEX/dWfHx8\n2Lt3b6Lsy/zLioRJdc5ePsvL815m8LrB7Dy+004xvQNjxoyhXLlyZM2alQIFCvDiiy8SERFx0216\n9uzJ8OHDPdp/Qta9lZsN47Ft2zZq165N7ty5yZUrFw8//LDXpwmtXr167GRCKZkVCZOqqCrP//Qc\nB35ZyuYD6wgeH8zdX91N51mdmfzbZI5fOO50xBRjwIAB9OzZkwEDBnDmzBnWrFnD/v37CQ4O5urV\nq/Fu4+QUnTf7Y6Bhw4bUrl2bv//+m3/++YeBAwfi5+eXhOn+K8VMZ+rJpBPJ/YZNOmRchm8YrmV7\nBej5+0urRkZqdHS0bj+2XQeuGagNJzZUv35+WvGb8tpjYQ9dvHexXoq85Gje5Pq7e+bMGc2WLZtO\nmzbtuufPnTunefLk0dGjR6uqap8+fbR58+batm1b9ff315EjR2qfPn20bdu2sduMHTtWAwMDNSAg\nQD/66KPrJh6Ku+6+fftURHTs2LFapEgRzZMnj/bt2zd2P+vWrdPKlStrjhw5tECBAvrSSy9pZGRk\n7HIR0T179vznvRw/flx9fHw0IiIi3vcaHh6uhQoV0k8++UQDAgK0aNGiOmHChNjlly9f1tdff12L\nFCmi+fLl0xdeeEEvXfr392bGjBn64IMPqp+fn5YoUUIXLFig7777rqZLl04zZ86s2bNn127dusVm\nHDJkiN5zzz1arFix2PccFRUVu7+goKDrJiGqUqWKvvrqq5ojRw4tXry4rlq1SseMGaOFCxfWvHnz\n6tixY2/yL3nnkw45/gWfGLfk+h/NJK1fjvyiAR/76fay+VSPHYt3nStHDumykhn1/dZ5tdI7d2n2\nDzJpnW+q6IAV/9Nfj/6q0dHRSZo5uf7uzp8/X319fa/78rqmQ4cO2qZNG1WN+ZLPkCGDzpo1S1VV\nL168qH369NF27dqpqurvv/+u2bJl01WrVmlkZKS+8cYbmiFDhuuKxLV1r31hPvfcc3r58mXdsmWL\nZsyYUXfs2KGqqhs3btS1a9dqdHS07t+/X++991796quvYnO5KxKqqiVLltQGDRrojBkz9O+//75u\nWXh4uKZPn17feOMNvXLlii5dulSzZs2qu3btUlXV7t27a+PGjfX06dN67tw5bdSokb7zzjuqqrp2\n7Vr19/ePfT+HDx/WnTt3qur1X/ZxM8adznTfvn3q4+Nz0yLh6+sbO4Xpe++9p0WKFNGXXnpJr1y5\nomFhYZo9e3Y9f/6823/LOy0SNnaTSRXOXj5Ly4lP8n/zlNIjfoSAgHjX881XgKq/RlD111/5cO1a\nTq1fxpKlK1i4vj9DHviGi5EXqVWsFsHFgqlVrBb5s+dP4ndyPfnglqMmeER7J6xf5vjx4wQEBODj\n898j0vnz52fTpk2xjytXrkzDhg0ByJQp03Xr/vDDDzRq1IjKlSsD8OGHHzJw4EC3rysi9OnThwwZ\nMlCuXDkeeOABtmzZQqlSpahQoULsekWKFOG5555j6dKlvPzyy7d8P0uWLOHTTz/ljTfe4M8//+Tx\nxx/n22+/jZ2ISET46KOP8PX1pVq1atSvX58pU6bw7rvvMmLECLZu3Yq/vz8APXr04KmnnqJv376M\nGjWKZ555hho1asR+NtcmaXLnnXfeid2XJ4oWLUr79u0BaNWqFZ988gm9e/fG19eX4OBgMmTIwB9/\n/EG5cuU83mdCWJEwKZ6q0mVmZ6ptjaBt60/gViNqZswIDz8MDz9MTl6iKdD08mXImJG9p/aycM9C\nZuycwSvzX6FQhgCCz+cluHR9qlXvSJbc+ZLkPV2T0C/3xBIQEMDx48eJjo7+T6E4cuQIAXGKcHzT\njl5z43SdmTNnJnfu3Dd9bXdThe7evZvXXnuNDRs2cPHiRa5evUrFihU9ej8FChSILU6HDh3i2Wef\npUOHDqxcuRKAnDlzXlfgAgMDOXz4MMeOHePChQvXvU50dHRs/8dff/1F/fr1PcpwjbvpSt2J+3lk\nzpwZ4LrPP3PmzLGfkTdYx7VJ8UZsGsFvJ7YxsOYA6Nr19nbiGoWzWM5idHmoCz+0/IFjbx5jRMUP\nyHk+mk9W9CPvF/mp+WI2Pu3+EBtnfE20Rifiu0heKleuTMaMGZk+ffp1z587d4558+ZRq1at2Odu\ndlbRtSlKr7l48SInTpy4rUwvvPACZcqUYc+ePZw+fZq+ffve1plrBQsWpGvXrvz222+xz506dYqL\nFy/GPj5w4AAFChQgICCALFmy8Pvvv3Py5ElOnjzJ6dOnY8/wKly4cOwsfDdy97nEfT5r1qwAXLhw\nIfa5o0ePJvg9eZMVCZOibTm6hXd/fpcpLaaSuW1HuMkXVkKl80nHo9VCea/fSpZ9GcHhnid4teHH\nHC6QnbY7PyVv/7y0ntaakZtGciDCNbfwhQtwG19cyY2fnx+9evWiW7duLFiwgKtXr7Jv3z5atWpF\nkSJFaNu2rUf7ad68OT/99BNr1qwhMjKSPn363HT9m33pnz17Fj8/P7JkycKOHTv45ptvPMpw+vRp\n+vTpw549e1BVjh8/zqhRo2IPgV173d69exMZGcny5cuZM2cOLVu2RER49tln6d69O8eOHQNiWiJh\nYWEAPPPMM4wePZolS5agqhw+fJidO3cCMS2AW123ERAQQMGCBfnuu++Ijo5m1KhRbouOJ5+RN1iR\nMCnW2ctnaTG1BV/W/pLSAaW9/nrZs+WiQd3uDHxrCdvfPsCm5zYRUjyERX8uouLwipQeXJpuHzzK\nrEo5OdOkLnz8MSxcCKdPez2bN7z55pt88sknvPHGG/j7+1O5cmUCAwNZtGgRvr6+Hu3j3nvvZdCg\nQbRq1YoCBQrg5+fHXXfd5Xb+hJtNFdq/f38mTJiAn58fXbp0oXXr1jfd9poMGTKwb98+goOD8ff3\np1y5cmTKlInRo0fHrpM/f35y5sxJgQIFaNeuHcOGDeOee+4B4LPPPqNEiRJUqlSJHDlyEBISwq5d\nuwB4+OGHGT16NN27d8ff35+goCAOHIj5g+GVV15h6tSp5M6dm+7du7vNOGLECD7//HMCAgLYvn07\nVapUuelnerPPyBu8Pp+EiNQBviSmII1U1c9uWN4GeNv18Czwoqr+6lq2D4gAooFIVX3EzWtoUldX\n4yxV5anpT5EtQzaGN0yci7HuRLRGs+XoFsL2LGDhttms/XsjD0bmJnh/eoJX/8PDQ2eRvkat/2yX\n1uaTOH/+PDly5OCPP/4gMDDQ6TgALF26lHbt2sV+uac2dzqfhFc7rkXEBxgM1AQOA+tFZKaq7oiz\n2l6gmqpGuArKcOBaz2M0EKSqp7yZ06Q8IzYO57cjW1jbZYPTUQDwER/K5y9P+fzlefvxHlyIvMCK\nAysI2xPG85XCOLC+BdWPVSe4WDAhxUMonqu405GTzOzZs6lZsybR0dG8/vrrlCtXLtkUCHNr3j67\n6RFgt6ruBxCRSUBjILZIqOqaOOuvAQrGeSzYITFzg1+O/sK7c15nxY7KZH4ps9Nx4pXFNwshxUMI\nKR4CwNFzR1m0dxEL9y7ko2UfkSl9JoKLBTucMmnMnDmTdu3aAfDQQw8xadIkhxOZhPDq4SYRaQbU\nVtXnXI/bAo+oarwnNovIG0DJOOvvBU4DUcBwVR3hZjs73JRGnLl8hoe+uo/eP53hqQlboUgRpyMl\nmKry+7HfWbhnIa899lqaOtxkkl6yPtyUECJSHegEPB7n6SqqekRE8gALRWS7qq6Ib/u4Z00EBQUR\nFBTkxbTGCarKc1PaEbTpBE+9PyNFFgiI+c9Z9q6ylL2rLK/xmtNxTBoRHh5OeHh4grfzdkuiEtBH\nVeu4Hvcg5lLwGzuvywE/AHVUNd7zv0SkN3BWVb+IZ5m1JNKAYWu/5uupb7Em+6tk7v2R03ESRVrr\nuDZJ705bEt4+3r8eKCEigSKSAWgNzIq7gogUIaZAtItbIEQki4hkc93PCoQAv2HSpM1HNvPeop5M\n2VuBzO9/4HQcY9IMrx5uUtUoEXkJCOPfU2C3i0iXmMU6HHgfyAV8LTEn/F471TUv8KOIqCvnBFUN\n82ZekzyduXyGltNaMrDRUEr1aAnxjCeUUgUGBnr9PHeTtt3pmWRev04iKdjhptRLVQn9IRT/jP4M\nazjM6TjJTsSlCJbsW0LYnjAW7l3ImctnYgcoDC4WTEG/grBhA9StC/PmwUMPOR3ZJBOeHm6yImGS\ntaEbhjJ0w1BWP7OazL7J83TX5GTf6X0s3LOQsL1h/Pznz+TLlo/GpRrzQUR5fLt1h9WrU2yHv0lc\nViRMirf5yGZCvgth5dMrKZm7pNNxUpyo6Cg2HdlEr/Be5Mqci/GnquNT7QlwDTdh0rbk0nFtzG05\nc/kMLUbXYWCJblYgblM6n3Q8XPBhpreczqEzh+iWbzPqmj/BGE9ZkTDJjqry7Ljm1Np8htAyLZ2O\nk+Jl9s3MrNBZrD20lveXvO90HJPCJJuL6Yy5ZuiyL9j521LWNBsJpb0/umta4JfRj/lt51N1dFVy\nZsrJ64+97nQkk0JYS8IkK5sPbaTXwneYmq41mVp7NmeB8UxAlgAWtlvIoHWDGLlpZMy8F7eYu8AY\na0mYZCPiUgQtvg1h0LZA7pkU7zBd5g4V8ivEwnYLeWLME/ifuUTzlh/AggVQvrzT0UwyZWc3mWRB\nVWk1rRW5j5/nm4bDIIHzAJuE2XJ0CyHfhTAu4Dlqvz8m5tRY+8zTFDsF1qQoX6//mhGbRrD6mdVk\nSp/p1huYO7bywEqenPwkM64257Epq2H5csiWzelYJolYkTApxqYjm6j9XW1WPb2Ke3LbOfxJacEf\nC2g/oz1he6vwwIErMHMmpEvndCyTBOw6CZMiRFyKoOXUlgyuO9gKhANql6jN4LqDqVd8LbuL54CT\nJ52OZJIZ67g2jlFVOv/QgZDiIbQq28rpOGlWi/taEHE5gpDlfVme8TLWM2HisiJhHPP1nD7sWT2H\n8e/udTpKmte5QmdOXTxFyPgQlnVaRkCWAKcjmWTCDjcZR2zcs4I+q/oypfT7ZMpf2Ok4Bnizyps0\nKd2EOt/V4czlM07HMcmEdVybJBdx8TQV+ham34nytPx6Kdh8CsmGqtJ1ble2HdvGvDZzyXzhCuTI\n4XQs4wV2dpNJllSVFp9VJO/v+xky9ABkzep0JHODaI2m7fS2nN2/m+kTo/Bdutz+nVIhO7vJJEtD\nVv0fe49uY8A74fbFk0z5iA9jm4xF8+WlU/UIotuEQlSU07GMQ6wlYZLMhsMbqDehHqs7Lqd4nlJO\nxzG3cDHyInXG1+b+NXsZlL0lMuALpyOZRGQtCZOsnL50mlbTWjGk3hArEClEZt/MzGrzE6sfDKDX\ngXEwdKjTkYwDrCVhvE5VaT61Ofmz5WdwvcFOxzEJdOz8MaoOr8Rz2zLz2tAtdkV2KmEtCZNsDF43\nmH2n9zEgZIDTUcxtyJM1DwufDmfgfecY9etYp+OYJGYX0xmv2rB5Dh8teIfVL/1CxvQZnY5jblNh\n/8KEtQsjaEwQfhn9aH5vc6cjmSRiRcJ4zekz/9ByUjO+zt6U4rmKOx3H3KGSuUsy96m51P6uNv4Z\n/QkuHux0JJMErE/CeIWq0uy9eyh4KopBg/eAjx3ZTC2uDTE+s+V0Khd+zP5tUyjrkzCOGjTsafaf\n/Yv+H6y2L5FUpkqRKox7chxNRtfh155POx3HeJn97zWJbv3KqXz851imNJ9Cxjz5nI5jvKBOiToM\nrPMldfmOP77p63Qc40XWJ2ES1elLp2m19g2+ue8tildr7HQc40WtKnXmzKmjBC/qxYrZxSjYINTp\nSMYLrE/CJBpVpdmUZhTyK8TAugOdjmOSyOfjn2fMhm9Z1iGcgAqPOx3HeMjTPglrSZhEM3DtQA5E\nHOD7Zt87HcUkobfaDeXUyUPUndqExfftxS+jn9ORTCKyloRJFOsOraPBxAas6byGYjmLOR3HJDFV\n5cU5L7DjxE7mPTWPTOkzOR3J3IKd3WSSzKljB2g1rRVDGwy1ApFGiQiD6w0hX7Z8tJrWisioSKcj\nmURiRcLcEY2KotPHD9MoshhNyzR1Oo5xUDqfdIxtMpbIqEienvU00RrtdCSTCKxImDvy1ccNOJT+\nAp+/NMvpKCYZyJAuA9NaTmP/6f28Mu9l7DBwymdFwty2dT8M5JOLYUx5bhEZM9kEQiZGFt8s/BT6\nEyuXT6T36PZOxzF3yIqEuS2ndm+l1apXGfZQH4qWetTpOCaZ8c/kz/ygb5n86/f838weTscxd8Dr\nRUJE6ojIDhHZJSJvx7O8jYhscd1WiEg5T7c1zlBVOo1pTOOclXmy+ftOxzHJ1F21m7Kw7Gd8uWIA\no5fZdTMplVeLhIj4AIOB2sB9QKiIlL5htb1ANVV9APgYGJ6AbY0DvlzzJYcDc/N5j8VORzHJXJHO\nrxOWsTPvzH+D6VsmOR3H3AZvtyQeAXar6n5VjQQmAdeN1aCqa1Q1wvVwDVDQ021N0lt7cC39VvRj\ncospZLD5IYwHSn04hLmHqvP8jM4s2rvI6TgmgTwqEiIyXUTqu/66T4iCwF9xHh/k3yIQn87AvNvc\n1njZyYsnaTWtFcMaDKNozqJOxzEphY8P5YfP4od2swn9IZQ1B9c4ncgkgKdf+l8DbYDdIvKpiCT6\nTPYiUh3oBFjfQzKkqnSa2YknSz/Jk2WedDqOSWkyZqRqsSDGNhlL40mN2fr3VqcTGQ95NHaTqi4C\nFomIPxDquv8XMAL4znU4KD6HgCJxHhdyPXcdV2f1cKCOqp5KyLbX9OnTJ/Z+UFAQQUFBN39TxnOq\n/N8nDTmS6zBTW0x1Oo1JwerdU4+BdQZSZ0IdlnZcSolcJZyOlGaEh4cTHh6e4O08HrtJRHIDbYF2\nwGFgAvA4cL+qBrnZJh2wE6gJHAHWAaGquj3OOkWAxUA7VV2TkG3jrGtjN3nRmv97jUbHBrG2268U\nzV/G6TgmFRi+cTifrviU5Z2WU9DPjiI7IVHHbhKRH4HlQBagoao2UtXJqtoNyOZuO1WNAl4CwoDf\ngUmqul1EuojIc67V3gdyAV+LyGYRWXezbT3JaxLPyWULaH3wK4aHDLICYRLNc2U70GW3PyFjanLi\nwgmn45ib8KglISL1VHXuDc9lVNXLXkuWANaS8A795x8a9wikxEMhfPHiTKfjmNSmRw/e/vs7llTO\nz+IOP5M9Y3anE6UpnrYkPC0Sm1S1wq2ec4oVCS+IimLAM2WYUuQcy3vvI0O6DE4nMqlNdDTaqiXP\nF9rC7gcLM/epuTbEeBJKlMNNIpJPRCoCmUWkvIhUcN2CiDn0ZFKp1QdX81nxo0x+eZkVCOMdPj7I\nuPF8vSond/35D62nteZq9FWnU5kb3LQlISIdgI7AQ8CGOIvOAmNUdbpX03nIWhKJ6+TFk5QfVp6B\ndQbSuLRdv2i87OhRrjz2KE1ezU9AoZKMaTIGnwRfkmUSKrEPNzVT1R8SJZkXWJFIPKpKo0mNKJmr\nJANqD3A6jkkrDh/mQm4/ak+sS/l85fmqzleI3PL7y9yBRCkSItJWVb8TkdeB/6yoql/cWczEYUUi\n8QxYNYCKby9MAAAb50lEQVSp26ayrJMdZjJJ7/Sl01QfW53GpRrTJ6iP03FSNU+LxK0uprs2SYDb\n01xN6rF65yI+X/U56zqvswJhHJEjUw7mPzWfqqOrkjNTTl6p9IrTkdI8jy+mS86sJXHnTowbSoXt\n3RnUcQqNSjVyOo5J4w5EHKDq6Kp8EPQBHR/s6HScVClRWhIictNB4FX15YQGM8lP9G9b6bDkZVrU\nCrUCYZKFIhczsCCiEdUX98Q/o7+NF+agWx1u2pgkKYxzzp7li/dqcuLRIvRr/a3TaYyJ4edH6dlr\nmNPoSerM7kL2jNmpVayW06nSJDvclJapsqpjTZ4MXM2613YQmCPQ6UTG/OvIEahUiWW9OtD85FBm\nhc6iUqFKTqdKNRLr7KYvVbW7iPxE/Gc3JYtjE1Ykbs+JpfMpP68RQ9p+T8OyzZyOY8x//for1KrF\nnBFv8vTu/ixqt4j7897vdKpUIbGKREVV3SgiT8S3XFWX3kHGRGNFIuGiNZpG3zeidM576F/3/5yO\nY4x78+ZBp058P603b67ry7JOyyiWs5jTqVK8ROm4VtWNrp9LRSQDUJqYFsVOVb2SKEmNIwasGsCJ\niyfo1+pHp6MYc3N168LcuYSWL09EJiF4fDDLOy2nQPYCTidLEzy94ro+MBTYAwhQFOiiqvNuumES\nsZZEwqw8sJKmU5qy/tn1FPEvcusNjElG+i3vx4StE1jacSm5s+R2Ok6KldjDcuwAGqjqH67HxYE5\nqlr6jpMmAisSnjt+4TgVhlVgSL0hNCzV0Ok4xiSYqvL2ordZun8pi9otsiHGb1OiTjoEnL1WIFz2\nEjPIn0lBoufPo8OohrS6r5UVCJNiiQif1fqMB/I+QJPJTbh09ZLTkVK1W3VcN3XdDQYCgSnE9Em0\nAA6o6oteT+gBa0l4YN8+Pu9Slhl1i7K02yZ80/k6nciY27dnD1Hbf6fNpQlcibrC1BZTSe9zq8u+\nTFyJ1ZJo6LplAv4GngCCgGNA5jvMaJLKpUusfLYOA6oIkzrNsQJhUr6zZ0n3dGfGF3iJS1cv8cys\nZ4jWaKdTpUp2MV0acPyFDlTINZWv20+mgR1mMqnFnDnw7LOcX7qI2sufo2L+inxZ50sbYtxDid1x\nnQl4BriPmFYFAKr69J2ETCxWJNyLHjeWBitepGyjznze4Cun4xiTuAYNgqFDOf3zXIJ+bMyTpZ+k\nd1Bvp1OlCIndcT0eyAfUBpYChbCO6xThf1HLiLi/JH3r9nc6ijGJr1s3qFmTHG07syB0LhO2TuCr\nNfbHUGLytCWxWVXLi8ivqlpORHyB5aqaLAZSsZZE/FYcWEHzKc1Z/+x6CvsXdjqOMd4RFRVz6KlR\nI/af3k/V0VX5qPpHdHiwg9PJkrXEmnTomkjXz9MiUhY4Ctx1u+GM9x07f4zQH0IZ2WikFQiTuqVL\nB41ihpELzBFIWLswqo+tjn8mf5qUbuJwuJTP0yIxXERyAu8Ds4iZqe59r6UydyRao2k/oz1tyrah\nfsn6TscxJkmVDijN7NDZ1J1Ql+wZslOzWE2nI6VodnZTanP+PJ9uHsRPu34ivEO4ne5q0qyl+5bS\nfGpzZofO5tFCjzodJ9lJ1I5rEcktIoNEZJOIbBSRL0XEBk1Jbv7+m+U1ivPlygFMajbJCoRJu37/\nnSfkbkY3Hk3jSY357Z/fnE6UYnl6dtMk4B+gGdAcOA5M9lYocxsiIznWvhlt6l5gVNOx1g9h0raf\nf4YGDWiQtypf1P6COt/VYe+pvU6nSpE8PbvpN1Ute8NzW1U1Wcz+keYPN0VHc7VDOxrkWciDdTvy\nafDnTicyxlmq0LUr7N0Ls2fzzeYR9F/dnxWdVpA/e36n0yULiX2dRJiItBYRH9etJbDgziKaxKJv\nvcmLmRYTVa4sH9Xo63QcY5wnAgMHxvzs1o0XHnqeZ8o/Q/D4YE5cOOF0uhTlVgP8nSVmQD8BsgLX\nBkfxAc6pqp/XE3ogTbckjh6ld89KzHkkB0ueXm7DJhsT15kz8Pjj0LEj+uqrvL3obZbsW8Li9ovx\ny5gsvr4ckygtCVXNrqp+rp8+qpredfNJLgUirfv6r+lMrODL3PZhViCMuZGfH8yeDXnzxg4x/lD+\nh2j4fUMuRF5wOl2K4PEpsCLSCKjmehiuqrO9liqB0mpLYtq2abwy/xWWd1puc/4a46FojabDjA4c\nv3CcGa1mkDF9RqcjOSKxB/j7FHgYmOB6KhTYoKo97yhlIkmLRWLJn0toNa0VYe3CeDDfg07HMSZF\nuRp9lZZTW+IjPkxqPilNzkWR2EXiV+BB1ZgB20UkHbBZVcvdcdJEkKaKxNmz/HJ+DyHjQ5jcfDLV\ni1Z3OpExKdLlq5dpNKkR+bLlY3Tj0fiIp+fxpA6JfXYTQI449/0THsncsYMH2Vu5NPXH1mZIvSFW\nIIy5XWvXkvGPP5necjp7T+3l5Xkvk2b+0EwgT4tEP2CziIwRkbHARsDOtUxKJ0/yT6Oa1G55mXdr\n9KbFfS2cTmRMyrV7N9SsSda9fzE7dDZrDq7hncXvOJ0qWbplkZCYaZ5WAJWA6cAPQGVV9eiKaxGp\nIyI7RGSXiLwdz/JSIrJKRC6JyGs3LNsnIltEZLOIrPPoHaVGFy5wtkk96jWIILTqi7z4cLKYWtyY\nlKttW+jXD2rWxP/Pw8xvO59Zu2bRb3k/p5MlO7fsrVFVFZG5rqurZyVk5yLiAwwGagKHgfUiMlNV\nd8RZ7QTQDYhvTN9oIEhVTyXkdVOVq1e50roFTSvvp8JDDfkg6AOnExmTOrRvH3OxXc2aBCxaxMJ2\nC6k2uhrZMmSj26PdnE6XbHjapb9JRB5W1fUJ3P8jwG5V3Q8gIpOAxkBskVDV48BxEWkQz/ZCwvpN\nUp3ovw7QsdR2sj74MF83+Mbm7zUmMbVrF1Mo6tShwPbtLGq/iGqjq5E9Y3Y6PtjR6XTJgqdF4lGg\nrYjsA84T8+WtHpzdVBD4K87jg8QUDk8psFBEooDhqjoiAdumeKrKazsH8te9BQlrPjlNnqZnjNe1\nbQvVqkHWrNxNVsLahVFjbA2y+ma1vj88LxK1vZrCvSqqekRE8hBTLLar6gqHsiS5z1d+zuI/F7Os\n4zIy+2Z2Oo4xqVeRIrF3SweUZt5T8wj5LoSsGbJS7556DgZz3k2LhIhkAp4HSgBbgZGqejUB+z8E\nFInzuJDrOY+o6hHXz2Mi8iMxrZB4i0SfPn1i7wcFBREUFJSAmMnP6M2j+WbDN6x8eiU5M+d0Oo4x\nacoD+R5gZuuZNPq+EVNaTCHo7iCnI92x8PBwwsPDE7zdrQb4m0zM/NbLgbrAflV9xeOdx1x0t5OY\njusjwDogVFW3x7Nub2IGDRzgepwF8FHVcyKSFQgDPlDVsHi2TT0X0/31F7MvbqHzrM6EdwyndEBp\npxMZkzZFRvLzweW0mtaKOW3m8EjBhBwpT/4S5YrruHNGiEh6YJ2qVkhgkDrAV8R0QI9U1U9FpAsx\nfRrDRSQvsAHITszZTOeAe4E8wI/E9EukByao6qduXiN1FIn161n9dAiNnvJhdtu5NuWiMU45dQoe\nfRSmTWN2pgM8M+sZFrZbSLm8yWKQiUSRWEViU9yicOPj5CJVFIldu9jWpArVn4pkTMvvqXtPXacT\nGZO2TZ4M3bvDggVM9tnOqwteJbxjOCVzl3Q6WaLwtEjcquP6ARE5c22fQGbX42tnN9lw4YnhyBH+\nalqLum2i+V+DgVYgjEkOWrWKOT22dm1aLVjAueofETw+mGUdlxGYI9DpdEnmpkVCVdMlVZA0KyKC\nkw1rUafFZbrV6En7B9o7ncgYc03LlrGF4pn58zlb6VVqja/F8k7LyZctn9PpkoSdeO+wC8cO07Dh\nOepWassbj73hdBxjzI1auK6V+PNPujfpztnLZwkeH0x4h3ByZ8ntbLYk4PGkQ8lZSu2TuBp9lScn\nP0mOTDkY22Rsmhuq2JiUSFV5e9HbhO8LZ1H7RSl2GtREnU8iuUuJRUJV6TyrM4fOHuKn0J/wTefr\ndCRjjIdUlRfnvMi249uY99Q8svhmcTpSgnljPgmTiN77+T22/rOVaS2nWYEwJoUREYbUH0IR/yI0\nm9KMK1FXnI7kNdaSSGobNzIwcgVDNnzNik4ryJM1j9OJjDG34+efueqXjRZ7+pFO0qW4aVCtJZEc\nTZvGpFeD+XzFZyxou8AKhDEp2ZkzpK/fkElF3+LM5TN0ntWZ6JgZnlMVKxJJJTycRZ905uW6MLfd\nfO7OcbfTiYwxd6JJExg2jIwNm/DjPe+z59SeVDkNqhWJpPDLL2x68UnaNBemtZmRqi7tNyZNa9IE\nhg8na6NmzC71IasPrubdn991OlWiSjkH0FKqvXv5I7Q2DdoKw54cSbXAak4nMsYkpsaNQQT/Vu1Z\nsGkFT0ytR/YM2elZtafTyRKFtSS87GjkKWq3g961+/FkmSedjmOM8YZGjWDzZgLyBLKw3UK+3fwt\ng9cNdjpVorCWhBeduXyGuuGd6VClK10e6uJ0HGOMNwUEAFAgewEWt18cO192Sp8G1U6B9ZLLVy9T\nb2I9SuUuxZB6Q2xuamPSmB3Hd1B9bHUG1R1E83ubOx3nP+yKawdFRUcR+kMo0RrN5OaTSedj4yQa\nkxb98udqav/QhNGNRye7aVDtOgknqKIzZ/LK/Jf55/w/fNf0OysQxqRVp07x4BOtmHHfR3SY0YGl\n+5Y6nei2WJFITP368cmELqzYt5yZrWeSKX0mpxMZY5ySMycMHUrlDu8xuUwvWkxtwbpD65xOlWBW\nJBLLt9/ybfgXjKyUkXntFuCfyd/pRMYYp9WrB2PHUqPTh4ws+QYNv2/I1r+3Op0qQaxIJIaZM5k5\n8k3er+nDgg6LyJ89v9OJjDHJRd26MH48DZ/tz1fFulJnQh12ndjldCqPWcf1nVq1ihXP16NpqA9z\nO4TxUIGHnMlhjEneFiwAYGSeg3y47EOWd1pOEf8ijsVJrDmuzS38ljOSZq19+K7lJCsQxhj3atcG\n4Bng7JWz1BpXi2WdliX7aVDtcNMd2H96P3Xnt+X/Gg4mpHiI03GMMSlE90rdaVeuHcHjgzl58aTT\ncW7KDjfdpuMXjlN1dFW6VOxC90rdk/S1jTEpn6ry1sK3WLp/qSPToNp1El50/sp5GkxsQONSja1A\nGGNui4jw+eVqVCQ/Db9vyIXIC05HipcViYSIjCRy8EBaTm1B6YDS9KvZz+lExpgUTLJkYcg7Kyl8\nMUOynQbVioSnoqPRpzvReef/EIQRDUfYeEzGmDtTsyY+kyYzpvdmMp06R5sf2nA1+qrTqa5jRcJT\nb79Nj3RL2HV/Aaa0nIpvOl+nExljUoOaNUk/aQqTPtrOmaP7k900qFYkPNG/P1/8MZ5ZFbIyu+1c\nsvhmcTqRMSY1qVGDjN9P4ceBf/PHsZ3JahpUO7vpVn78kQkDn6Vn/YyseHa1oxe/GGNSuXPniEgf\nRY1xNahdvDaf1PzEay9lZzclkrBiymt1YF77MCsQxhjvypYN/0z+LGi7gJk7Z9JvufMnx1iRuIn1\nh9bTNux5pofO5L677nM6jjEmjQjIEpBspkG1YTnc2HViF40mNeLbRt9SpUgVp+MYY9KYAtkLsKjd\nIqqNetzRaVCtSMTj8NnD1P6uNh9X/5hGpRo5HccYk0YVlZwsHBlJ9cuvky1DNkemQbUiEdepU0T0\n/5i6hRbxbIVneabCM04nMsakZTlyUHr4dOY914gQ7UwW3yxJPg2q9Ulcc/EilxrXp7HPFJ64uxo9\nH+/pdCJjjIHHH+fBET8xc5LQYUqbJJ8G1U6BBbh6lahmTWlZ8hfSP1KJ75tPwkesfhpjkpGVK1nc\nrT6hzWB2xzAeKfjIHe0u2ZwCKyJ1RGSHiOwSkbfjWV5KRFaJyCUReS0h2yYKVbTLc3QtsJmIsiUY\n9+R4KxDGmOSnShVqDp7LyHt7Juk0qF5tSYiID7ALqAkcBtYDrVV1R5x1AoBAoAlwSlW/8HTbOPu4\n/ZbEN9/wwYq+zHw8gPCnlyX5cL3GGJNQk36bxOthr7OkwxJK5i55W/tILjPTPQLsVtX9rlCTgMZA\n7Be9qh4HjotIg4RumxiGlrvC+MsZWNlugRUIY0yK0Lpsa85dOUfw+GCvT4Pq7eMqBYG/4jw+6HrO\n29t6ZPr26Xy45jMWtF9I3mx5E3PXxhjjVZ0rdObVSq9Sa1wtjp476rXXSTWnwPbp0yf2flBQEEFB\nQTddf+m+pTw/+3kWtF1A8VzFvRvOGGO8oHul7pzZuIrgryuxtNsmcmXO5Xbd8PBwwsPDE/wa3u6T\nqAT0UdU6rsc9AFXVz+JZtzdwNk6fREK2TVCfxJajWwgeH8yk5pOoUbTG7bw1Y4xJFnTNGt76tAZL\nqxRi0csbPD5snlzObloPlBCRQBHJALQGZt1k/biBE7rtre3ezZ/Pt6b+xPoMrjfYCoQxJsWTSpX4\nvOcSKqz7i4aDH0v0aVC9fp2EiNQBviKmII1U1U9FpAsxrYLhIpIX2ABkB6KBc8C9qnouvm3dvMat\nWxJHjnCsRiWqPHWJbsHv0e3Rbon1Fo0xxnHRa9fQ/qsgTj5clhkvryJDugw3Xd/TlkTauJguIoJz\n1atQo/FpQqp25OMaHyddOGOMSSKRa1fTYkQw6WvUYlLraaT3cd/tbEXimkuXuFI3hIaP7aPwo8GM\naPitzU1tjEm1Ll84S8MfmlEgewFGNR7l9uLg5NIn4bjozz/j6Qf3k6nsgwxtMMwKhDEmVcuYJTs/\ntvqRP07+wSvzXrnjaVBTdZFQVd6oeIJ9ZQsyqfnkmza9jDEmtciaIStz2sxh1cFVvPvzu3e0r1Rd\nJPqv6k/Y/p+Z1WY2mX0zOx3HGGOSjH8mf+Y/NZ8ZO2bQb06P295Pqv3TetyWcQxeP5iVT6+86QUm\nxhiTWuXJmodFTaZTdUBZsh86zkvPfZvgfaTKlsTc3XN5a+FbzH9qPoX8CjkdxxhjHFOgYGkWNZ3B\nZ7tHM/bbhJ/6n7paEuHhrBn3CR3L/MKs0FmUyVPG6UTGGOO4oo83YKFOo/pPzck6KiPNn+7v8bap\np0j88gvbuzSlSQcY0+Q7KhWq5HQiY4xJNkpXfZK5UROoPb8NWcdm8ni7VFMkDrWoQ91O6fms3v+S\nfA5YY4xJCcoHtWamXqHRmlc83ibVXExX9qMCtKv+Cm9VecvpOMYYk6wt3ruYWsVrpa0rrl+d/yoD\nQgbYxXLGGOOBNDcsR1R0lM1NbYwxHkpzw3JYgTDGmMRn36zGGGPcsiJhjDHGLSsSxhhj3LIiYYwx\nxi0rEsYYY9yyImGMMcYtKxLGGGPcsiJhjDHGLSsSxhhj3LIiYYwxxi0rEsYYY9yyImGMMcYtKxLG\nGGPcsiJhjDHGLSsSxhhj3LIiYYwxxi0rEsYYY9yyImGMMcYtKxLGGGPcsiJhjDHGLSsSxhhj3LIi\nYYwxxi2vFwkRqSMiO0Rkl4i87WadgSKyW0R+EZHycZ7fJyJbRGSziKzzdlZjjDHX82qREBEfYDBQ\nG7gPCBWR0jesUxcorqr3AF2Ab+IsjgaCVLW8qj7izaypRXh4uNMRkgX7HP5ln8W/7LNIOG+3JB4B\ndqvqflWNBCYBjW9YpzEwDkBV1wL+IpLXtUySIGOqYv8JYtjn8C/7LP5ln0XCefsLuCDwV5zHB13P\n3WydQ3HWUWChiKwXkWe9ltIYY0y80jsd4BaqqOoREclDTLHYrqornA5ljDFphaiq93YuUgnoo6p1\nXI97AKqqn8VZZyiwRFUnux7vAJ5Q1b9v2Fdv4KyqfhHP63jvTRhjTCqlqnKrdbzdklgPlBCRQOAI\n0BoIvWGdWUBXYLKrqJxW1b9FJAvgo6rnRCQrEAJ8EN+LePJGjTHGJJxXi4SqRonIS0AYMf0fI1V1\nu4h0iVmsw1V1rojUE5E/gPNAJ9fmeYEfXa2E9MAEVQ3zZl5jjDHX8+rhJmOMMSlbij691JML9dIC\nERkpIn+LyK9OZ3GaiBQSkZ9F5HcR2SoiLzudySkiklFE1rouRt3q6tdL00TER0Q2icgsp7M4KSEX\nKqfYloTrQr1dQE3gMDH9H61VdYejwRwgIo8D54BxqlrO6TxOEpF8QD5V/UVEsgEbgcZp8fcCQESy\nqOoFEUkHrAReVtU0O3qBiLwKVAT8VLWR03mcIiJ7gYqqeupW66bkloQnF+qlCa7Tgm/5j50WqOpR\nVf3Fdf8csJ3/XpuTZqjqBdfdjMT07aXMvwoTgYgUAuoB3zqdJRnw+ELllFwkPLlQz6RhInI38CCw\n1tkkznEdXtkMHAUWqup6pzM56P+AN0nDhTIOjy9UTslFwhi3XIeapgGvuFoUaZKqRqtqeaAQ8KiI\n3Ot0JieISH3gb1crU1y3tKyKqlYgpmXV1XXIOl4puUgcAorEeVzI9ZxJ40QkPTEFYryqznQ6T3Kg\nqmeAJUAdp7M4pArQyHUs/nuguoiMcziTY1T1iOvnMeBHYg7fxyslF4nYC/VEJAMxF+ql5TMW7K+j\nf40CtqnqV04HcZKIBIiIv+t+ZiAYSJMd+Kr6jqoWUdVixHxX/Kyq7Z3O5QQRyeJqaRPnQuXf3K2f\nYouEqkYB1y7U+x2YpKrbnU3lDBGZCKwCSorIARHpdKttUisRqQI8BdRwnd63SUTS6l/P+YElIvIL\nMf0yC1R1rsOZjPPyAitcfVVrgJ9udqFyij0F1hhjjPel2JaEMcYY77MiYYwxxi0rEsYYY9yyImGM\nMcYtKxLGGGPcsiJhjDHGLSsSxriISJTruopr11e8lQSvOVxESnv7dYy5XXadhDEuInJGVf0SeZ/p\nXBd+GpMiWUvCmH/FO6yJiPwpIn1EZKNropaSruezuCZ8WuNa1tD1fAcRmSkii4FFEuNrEdkmIgtE\nZI6INHWtu0REKrjuB4vIKhHZICKTXfO8IyKfishvIvKLiHyeJJ+EMS5WJIz5V+YbDje1iLPsH1Wt\nCAwF3nA99y6wWFUrATWA/q4xkgDKA01VtTrQFCiiqvcC7YHKN76wiOQG3gNqqupDxEyW9JqI5AKa\nqGpZVX0Q+DjR37UxN5He6QDGJCMXXMMnx+dH18+NwJOu+yFAQxF50/U4A/+OTLxQVSNc9x8HpgKo\n6t8isiSe/VcC7gVWiogAvsSMxxUBXBSRb4E5wOzbemfG3CYrEsZ45rLrZxT//r8RoJmq7o67oohU\nAs4ncP8ChKnqU/9ZIPIIMdP0tiBmUMuaCdy3MbfNDjcZ86+EDrW+AHg5dmORB92stxJo5uqbyAsE\nxbPOGqCKiBR37SuLiNzjGso5h6rOB14D0vQc5ibpWUvCmH9lEpFNxBQLBear6ju4n+7yI+BLEfmV\nmD+49gKN4lnvB2L6LH4nZsrdjcQcRuLavlX1uIh0BL4XkYyu598DzgIzRSSTa/1X7+gdGpNAdgqs\nMUlARLKq6nlXR/RaYqaP/MfpXMbcirUkjEkas0UkBzEd0h9agTAphbUkjDHGuGUd18YYY9yyImGM\nMcYtKxLGGGPcsiJhjDHGLSsSxhhj3LIiYYwxxq3/BwqGCBVMuSBSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xba38d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot1, = plt.plot(lc/float(sum(lc)), 'r--', label='Assigned energies')\n",
    "plot2, = plt.plot(prob,'g',label='Original Spectrum')\n",
    "plt.xlabel('Energies')\n",
    "plt.ylabel('Probability')\n",
    "plt.legend(handles=[plot1,plot2])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}